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Abstract. The framework of third quantization - canonical quantization in the 

Liouville space - is developed for open many-body bosonic systems. We show how to 

rl^ ' diagonalize the quantum Liouvillean for an arbitrary quadratic n— boson Hamiltonian 

with arbitrary linear Lindblad couplings to the baths and, as an example, explicitly 
work out a general case of a single boson. 
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1. Introduction 

Liouville spaces of operators have proven very useful in the analysis of many-body 
problems p] , and in particular in optical spectroscopy [2] . In such approaches one equips 
the Liouville space of quantum operators (representing observables or density matrices) 
with an inner product such that the Liouville space in fact becomes a Hilbert space. 
Therefore, for many body problems, a second quantization over such spaces is tempting 
and indeed was presented by one of the authors [3], and coined as third quantization, 
for the finite dimensional case of Fermi operators. It was successfully used to explicitly 
and efficiently solve situations with quadratic (or quasi-free) Hamiltonians and linear 
coupling to an environment via the Lindblad operators [3j H] or via the Redfield model 
[5] . In particular, the method already provided exciting new physics, such as a discovery 
of a quantum phase transition far from equilibrium in XY spin chain p[] . 

The present communication deals with developing a similar framework for 
bosons. This is of considerable importance for several reasons. First, there exists a 
somewhat dated but vast literature on harmonic oscillators and more general quadratic 
Hamiltonians reported in two books [HI [7] essential to the development of quantum optics 
and many other fields. Second, oscillators and linear couplings to the environment are 
widely used in the theory of decoherence [8|. Finally, the development of technique 
of third quantization over infinite-dimensional boson spaces seems essential for a future 
development of a supersymmetric theory of open many-body systems involving fermions, 
bosons, their interaction, and coupling to the environment. For quadratic interactions 
and linear coupling to the environment closed form solutions will likely emerge, whereas 
for non-linearly interacting systems perturbative and non-perturbative methods of 
many-body physics can be used in conjunction with our approach. Though many 
aspects of the open bosonic problems with linear environment coupling operators have 
been treated before, the elegance and versatility of the third quantization treatment 
distinguishes this treatment from others. 

2. Preliminaries 

Let us consider a Hilbert-Fock space % of n bosons. The representation of "H can 
be generated by a special element ipo £ %, called a vacuum pure state, and a set 
of n unbounded operators over H, a±,...,a n , and their Hermitian adjoints a\, . . . , a) n , 
satisfying canonical commutation relations (CCR) 

[%■ , 4] = S J,k , [a,- , a*] = [a] , 4] = 0. ( 1 ) 

Let us define a pair of vector spaces /C and /C', such that /C contains trace class operators, 
for example density matrices, and K! contains unbounded operators which we need as 
physical observables. After choosing a specific space of observables /C, we define the 
space /C as a subspace of trace class operators over %, such that p G /C if and only if Ap 
is trace class for any A e KJ '. In this sense K! is a dual space to K,. 
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For instance, we may chose /C as a linear space of all (unbounded) operators whose 
Weyl symbol, i.e the phase space representation of the operator, is an entire function on 
C 2n . Then /C must be limited to operators with finite support in the number operator 
basis, i.e. to operators which have a finite number of non-vanishing matrix elements 
in this basis. Such a constraint on density matrices may be too restrictive for certain 
applications. We shall show later [using eq. ©] that it can be relaxed by restricting KJ . 

Since the development in the paper will be purely algebraic, we shall adopt Dirac 
notation and write an element of K! as ket \p) and an element of K, as bra (A\, such that 
their contraction gives the expectation value of A for a state p, 

(A\p)=trAp. (2) 

We keep distinct styles of brackets to emphasize the manifest difference between the 
spaces from which the ket and the bra have to be chosen. 

If b is any of a,j, a], then for each p G /C and A e K! , bp, pb and Ab, bA are also 
elements of /C and K.' respectively, so one can define the left and right multiplication 
maps b and b over K, 

b L \ P ) = \bp), b R \p) = \pb). (3) 

The action of their adjoint maps on KJ follows from definition (T5]) and cyclicity of trace, 

[A\b h = {Ab\, (A\b R = (bA\. (4) 

With a slight abuse of notation, we can also write (b ) = b and (b ) = b . 
Let us now define the set of 4n maps a v j, a' u j,j = I, ■■■ ,n,u = 0,1, 



QqJ — 0,j ; CLq j — Q'j — O' j 



a\ j = a) j , a[j = af — a}j. (5) 

satisfying the unique properties: (i) almost-CCR 

fo>j> a kfc] = 5 w 8 J,k> \Kv a»,k\ = K,j; o^] = 0, (6) 

(ii) a' V j left-annihilate the identity operator 

(l|aV = (7) 

and (iii) a V j right-annihilate the vacuum pure state \po) = l^oX^ol 

fli/j|po) = 0- ( 8 ) 

Writing a In component multi-index m = [m u j e Z + ; v 6 {0, 1}, j 6 {1 . . . n}) T we can 
define convenient dual-Fock bases of the spaces /C, K! as 

(fit \m Ulj ,f \ m ■ 



"■J 



vj- 



whose bi-orthonormality (m'|m) = 5rn',m is simply guaranteed by almost-CCR (jSJ). Here 
and below, x = (xx,x 2 , ■ ■ -) T designates a vector (column) of any, scalar-, operator- or 
map- valued symbols, x-y = Xiyi + x 2 y2 + ■ ■ ■ designates a dot product, and bold upright 
letters A shall be used for complex valued matrices. 
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The explicit construction of the bases ([2]) allows us to enlarge the space K and 
restrict the space KJ as compared to the above example. Namley, we now identify the 
space K with the I 2 Hilbert space of vectors of coefficients {am}, K,3 \a) = J2 m a m\ui), 
and the space K! with the I 2 Hilbert space of vectors of coefficients {Sm}, K? 3 (S\ = 
Em^mH- Then, clearly by Cauchy-Schwarts inequality, |trSVj| = | ^2 m Smp~m\ < °° 
and hence /C and K! are dual in the required sense. 

3. Explicit solution of the Lindblad equation for quadratic bosonic systems 

Our goal is to present an exact solution of the master equation of an open n-particle 
system, say of the Linblad form 

^ = tp := -i[H, p] + J2 ( 2L ^ L l - i L i L ^ p}) ( lo ) 

where H is a Hermitian operator (Hamiltonian), {x,y} := xy + yx, and L M are arbitrary 
(Lindblad) operators representing couplings to different baths. We are going to describe 
a general method of explicit solution of f llOp for an arbitrary quadratic system of n bosons 
with linear bath operators 

H = a) ■ Ha + a ■ Ka + a f • Ka f (11) 

L» = l_,-a + k ll -(J. (12) 

Note that the hermiticity of the Hamiltonian and CCR ([T]) implies that the matrix H is 
Hermitean H^ = H r = H and the matrix K is symmetric K = K T . For bound oscillator 
systems, one can always choose K = 0, however, in order to be able to describe freely 
moving or inverted oscillator modes, we shall keep the general form. 

Following the rules of the algebra ( I3f4"f5]) we straightforwardly express the 
Liouvillean ffTUl) in terms of the canonical maps 



L •"• , R A , L " L •*> R * , R 

x - ■ '" 1 r T T T T\ 



C= ~iH L + iH R + Y, 2L„ 4 - 4 ^ - L» L\ (13) 

= ~ i^o ' Hoq + i«i ■ HS-l + ia^ ■ K(2a + oQ — ia' ■ 0.(2(1! + a' ) 
+ a ■ (N-M)oo + a x ■ (N-M)oj + a ■ (L T -L)o 1 + aj • (L T -L)a 
- «o ' La' - ai • Lai + 2a' • Nai (14) 

where we define n x n matrices 

M := ^ ® ^ = Mt, N := £ A; M ® ^ = N+, L := ^ ® V ( 15 ) 

M M /^ 

The fact that the primed maps a' always appear on the left side in each term manifestly 
expresses the trace preservation of the flow (TTOl) . i.e. (1|£ = 0. 

Note that adding inhomogeneous terms (linear forces) to the Hamiltonian H — > 
H + f ■ a + / ■ a) and the Lindblad operators L M — > L M + A M 1, where / is a c- vector and 
A^ are c-numbers, simply results in an extra linear term in the Liouvillean £ — > C + g-a, 
which must be removed by a canonical transformation affecting the annihilation maps 
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a u j only, namely a — > a + si, a — > a and thus preserving the important relations (16"]T7|) . 
while the right vacuum state \po) (which is unimportant for the following discussion) is 
trivially shifted. 

Writing An vector of canonical maps as b = (a, a) T = (a , ai,a$,a.i) T the Liouvillean 
(114p can be compactly rewritten in terms of a symmetric form 

C = b-Sb-S l (16) 

where S is a complex symmetric An x An matrix which can be conveniently written in 
terms of two In x In matrices 



namely 



and 



-X 

-X T Y 



1 / iH - N + M -2iK - L + L r 



(17) 



X : ~ 2 V 2iK - L + 17 -iH - N + M ' (18) 



Y: — ( o. T r w T T r)=Y T (19) 



1 /-2iK-L-L T 2N 

2 V 2N T 2iK - L - L J 

The scalar So ( fT6l) stemming from reordering of maps equals So = tr X = tr M — tr N. 
For the rest of our discussion we shall assume that the matrix X is diagonalizable 
[9], i.e. it is similar to a diagonal matrix 

X = PAP\ A = diag{A,...,/3 2n } (20) 

where the In complex eigenvalues /3j shall be named as rapidities in analogy to the 
fermionic case [3]. Let J = ia y <8> l 2n denote the symplectic unit, satisfying J 2 = — l 4n . 
It is straightforward to check that the matrix JS belongs to the symplectic algebra, 
namely J(JS) = — (JS) T J, and can be diagonalized as 

>T „ T-1-I1 I Izn — Z 



JS = V- 1 [(-A)©A]V, V=[P T ©P- 1 ] '" fl (21) 

V U i2n / 

with the eigenvector matrix V belonging to the complex symplectic group 

V T JV = J. (22) 

The 2n x 2n complex symmetric matrix Z = Z T is a solution of the continuous Lyapunov 
equation [TO] 

X T Z + ZX = Y. (23) 

It is known that the solution of (|23|) exists and is unique if no pairs of eigenvalues, 
rapidities fy, (3j> exist, such that fy+fiji = 0. We note that X, and also Y, are (unitarily) 
similar to real matrices, namely if U := -75(12 + i^x) ® ^2n then UXU -1 , and UYU -1 , 
are real matricesjj Thus the rapidities should come in complex conjugate pairs /3j,f3j. 

I In fact, if we were working in terms of Hcrmitian "coordinate" and "momentum" operators, 
aj + a!:,i(aj — a A, from the outset, then a very similar formalism could have been developed with 
matrices X and Y being automatically real. 
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Existence and uniqueness of the solution of ( 1251) and further, the decomposition ( )2"T|) . 
is thus guaranteed if all the rapidities lie away from the imaginary axis, Re/3j ^ 0. 
Diagonalizing X ( 1201) is then also the essential part of efficiently computing the solution 
Zof (E3D [n]. 

Defining the 4n normal master-mode maps (NMM) as (£, £ ) T := V6, or 

( = p T (a - Za ), f = P -1 ^' (24) 

which, due to symplecticity ( l2"2"j) . satisfy almost-CCR 

[Cr,£]=^ jS , [Cr,6] = [C,C] =0, (25) 

brings the Liouvillean to the normal-form £ = Jb-JSb — Sol = J( ■ [(— A) © A]£ — Sol 
and finally 

In 

£=-2j2PrC r Cr- (26) 

r=l 

In contrast to the (finite dimensional) fermionic case [3] , existence of a stable fixed 
point of the Liouvillean dynamics is here not guaranteed. Since our (central) system 
is infinitely dimensional, it can absorb excitations from the environment and amplify 
indefinitely. Indeed, this is signaled by at least one of the rapidities being to the left of 
the imaginary axis, 3j, Re/3j < 0. 

However, if we assume that all the rapidities lie to the right of the imaginary axis, 
Vj, Re f3j > 0, then we can make the following statements: 

(i) A unique non-equlibrium steady state (NESS) exists |NESS) G /C namely the "right 
vacuum state" of the Liouvillean 

£|NESS) = 0. (27) 

All the physical properties of NESS are essentially determined by the NMM 
annihilation relations 

(1|C = 0, C r |NESS) = 0. (28) 

(ii) The complete (point) spectrum of the decay modes of the Liouvillean is given in 
terms of a 2n component multi-index of super-quantum numbers m G Zi n 



A M = -2^m r /3 r , £16™) = A^je^), (6J£ = \rn{®m\ (29) 



where 



|e ZZL ) = n i 7=|NESS), (e„| = (11X1^4=. (30) 



(iii) The 2-point correlator of NESS is given by the solution of the Lyapunov equation 
( 1231) . If b — {a,a}) T designates a In column of canonical operators then 

tr :b r b s : p NESS = (:b r b s : |NESS) = (l|a r a s |NESS) = Z r .. (31) 
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Proof. Using purely algebraic manipulations (i) and (ii) can be proved in analogy to the 
fermionic case j3] and are in fact equivalent to the standard construction of the Fock 
ground state and quasi-particle excitations with a restriction of non-normality, (iii) The 
first two equality signs are just the definitions ( 12|5|) whereas the last equality follows 
after expressing a r in terms of NMM maps (I24J1 , and using (1281) . 

We note that NESS is a essentially a Gaussian state and Wick theorem can be used 
to express any higher-order correlator in terms of 2-point contractions ( )3~Tj) . 

4. Example: a general Linbldad equation for a single quantum oscillator 

Take the simplest case, namely a single quantum harmonic oscillator n — 1. We 
choose H = u;,K = 0,M = u = J2„ K? > 0, N = v = ^ \k^\ 2 > 0, 
L = w = J2u Ifjfcn e ^-- By construction the reservoir parameters should satisfy the 
positivity condition uv > \w\ 2 . The matrices ( TT8lfT9l) now read 

x = I( iw + M ° V Y=(-* V V (32) 

The rapidity spectrum (3± = ^(u — v ± ico) indicates that the problem has a stable 
fixed point if and only if u > v . Then, NESS is approached at an exponential rate 
~ exp[— {u — v)t] and the complete corelator (I3T|) . as given by the straightforward 
solution of the Lyapunov equation (|23|) (now a linear system of 3 equations), read 



//' 



tra p NE ss = Z X1 = — — = tr (at) 2 p NESS , 

u — V + 100 



tr a^pNEss = Z 12 = — - — • (33) 

u — v 

It is perhaps worth noting that this is not a Gibbs thermal state for any combination of 
Lindblad parameters, the fact which can be viewed as a manifestation of non-ergodicity 
of the harmonic oscillator. Furhtermore, since NESS is a generalized Gaussian state, 
expectation values of all higher moments of the field are calculated trivially. 

5. Discussion 

In the present communication we outlined a simple and mathematically consistent 
framework for diagonalizing quantum Liouvilleans for arbitrary bi-linear systems of 
bosons. Our method in some sense resembles the technique known as non- equilibrium 
thermo-field dynamics [T2| H5] , which has been recently applied to a similar problem 
of diagonalizing the Lindblad equation [H]. However, there is an important difference, 
namely in thermo-field dynamics, one treats a density matrix as an element of a tensor 
product Hilbert space, whereas the operator (observable) of choice is fixed in a given 
calculation. In our approach, we consider (dual) spaces of both simultaneously, so the 
canonical structure of our formalism is simpler and more transparent. This, in turn, 
is reflected in the ease of getting simple general results as outlined in the example in 
section H]) . 
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Notice that our formalism is flexible with respect to taking different realizations of 
the dual spaces K and /C, namely for some physics applications it may be advantageous 
to take larger space of observables K' at the expense of a smaller dual space /C, or vice 
versa. In analogy to the fermionic case [5] this formalism can be extended to explicitly 
time dependent problems. For example, for periodic time dependent problems [15], the 
Floquet representation of the correlation matrix leads to a discrete Lyapunov equation 
rather than a continuous one that appears here. 

The approach presented here differs essentially from the framework proposed earlier 
for fermionic systems [3] , where the density operators and the observables belong to the 
same Hilbert space. For bosonic systems, due to infinite-dimensionality, this is not 
possible and the symmetry between density operator space and observable operator 
space is broken. However, it may be desirable for certain purposes to have formally 
similar construction for bosonic and fermionic systems. Thus we point out an equivalent 
version of fermionic third quantization, which follows exactly the steps of the present 
communication, but starting instead from a set of fermionic operators Cj,Cp obeying 
canonical anti- commutation relations (CAR), 

{cj, c{} = S j>k , {cj, c k } = {c], c{} = 0, (34) 

introducing the dual spaces of density operators and observables, stating ff2l[3l4|) and 
defining the canonical adjoint fermionic maps 



c 0yj = c^, c' 0j = ct, -ct,"P, 



,R 



3 3 



cij = c^V, c' hj =cfV-cj, (35) 

satisfying almost-CAR 

{c^-, c^J = 5 Vtfl 8 jtk , {cyj, c M)fc } = {c'^j, c' M) J = 0, (36) 

and the properties (!7f8|) . The parity superoperator V is uniquely defined by the dual 
vacuum states, {\\P = (l\,V\po) = \po), and anticommutes with all elements of the 
adjoint- algebra {V,c} = {V,c} = 0. The difference to the more symmetric approach 
[3] is that now the canonical conjugate adjoint maps are not the hermitian adjoint maps 
c' u j 7^ cj, -, which is however of no consequence as we are anyway dealing with problems 
in which non-normal operators enter in an essential way. 
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